;; Tests of the naive B-spline functions
;; Copyright (C) 2023  Rocks Mazama
;; 
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;; 
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;; 
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <https://www.gnu.org/licenses/>.

(ns mazama.compgeom.b-spline-test
  (:require
   [clojure.test :refer [deftest is testing]]
   [mazama.compgeom.b-spline]))

(def EPSILON 1/10000000)
 
;;; Reference: https://floating-point-gui.de/errors/comparison/
(defn nearly== [a b epsilon]
  (or (== a b)
      (if (and (rational? a) (rational? b))
        false
        (let [absA (abs a)
              absB (abs b)
              diff (abs (- a b))]
          (if (or (== a 0) (== b 0) (< (+ absA absB) Float/MIN_NORMAL))
            (< diff (* epsilon Float/MIN_NORMAL))
            (< (/ diff (min (+ absA absB) Float/MAX_VALUE)) epsilon))))))

(defn- p== [p q epsilon]
  (cond
    (and (nil? p) (nil? q)) true
    (and (number? p) (number? q)) (nearly== p q epsilon)
    (and (coll? p) (coll? q)) (and (== (count p) (count q)) (every? identity (map #(p== %1 %2 epsilon) p q)))
    :else false))

(defn- v+
  "Add two vectors."
  [v w]
  (map + v w))

(defn- v*
  "Multiply a scalar and a vector."
  [k v]
  (map #(* k %) v))

(defn- check-domain [f knots]
  (is (every? identity (map f knots)))
  (let [delta 1/10000000
        k0 (first knots)
        kn (last knots)]
    (is (nil? (f (- k0 delta))))
    (is (nil? (f (- k0 (/ delta)))))
    (is (nil? (f (+ kn delta))))
    (is (nil? (f (+ kn (/ delta)))))))

(defn- check-spline
  ([f knots domain-range-map]
   (check-spline f knots domain-range-map EPSILON))
  ([f knots domain-range-map epsilon]
   (check-domain f knots)
   (doseq [[t v] domain-range-map] (is (p== (f t) v epsilon)))))

(defn- lerp
  "Parametric linear interpolation from `p0` to `p1`.
  `f+` is a function that adds two points together;
  `f*` is a function that multiplies a scalar and a point."
  [f+ f* p0 p1 t]
  (f+ (f* (- 1 t) p0) (f* t p1)))

(defn- shene-de-Boor
  "A helper function used in the example at the bottom of https://pages.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/de-Boor.html"
  [f+ f* control-points]
  (let [p1_0 (nth control-points 1)
        p2_0 (nth control-points 2)
        p3_0 (nth control-points 3)
        p4_0 (nth control-points 4)
        p4_1 (lerp f+ f* p3_0 p4_0 1/5)
        p3_1 (lerp f+ f* p2_0 p3_0 8/15)
        p2_1 (lerp f+ f* p1_0 p2_0 4/5)
        p4_2 (lerp f+ f* p3_1 p4_1 3/10)
        p3_2 (lerp f+ f* p2_1 p3_1 4/5)]
    (lerp f+ f* p3_2 p4_2 3/5)))

(deftest error-test
  (testing "no control points"
    (is (thrown? Error (mazama.compgeom.b-spline/b-spline + * [] [0 0 0 0 0.5 1 1 1 1]))))
  (testing "too few knots for the given control points"
    (is (thrown? Error (mazama.compgeom.b-spline/b-spline + * [2 4 6] [0 1 2])))))

(deftest degree-0-test
  (testing "{:knots 2, :points 1, :degree 0}"
    (let [knots [0 1]
          b-spline (mazama.compgeom.b-spline/b-spline + *
                                                      [-37]
                                                      knots)]
      (check-spline b-spline
                    knots
                    {0 -37
                     1 -37})))
  (testing "{:knots 3, :points 2, :degree 0}"
    (let [knots [-1 1 3]
          b-spline (mazama.compgeom.b-spline/b-spline + *
                                                      [73 -10]
                                                      knots)]
      (check-spline b-spline
                    knots
                    {-1 73
                     0 73
                     1 -10
                     3 -10})))
  (testing "{:knots 4, :points 3, :degree 0}"
    (let [knots [1 11 21 31]
          b-spline (mazama.compgeom.b-spline/b-spline + *
                                                      [10 -73 66]
                                                      knots)]
      (check-spline b-spline
                    knots
                    {1 10
                     11/2 10
                     11 -73
                     2099/100 -73
                     21 66
                     31 66}))))

(deftest degree-3-test
  (testing "http://nurbscalculator.in/ 1"
    (let [knots [0 0 0 0 1 1 1 1]
          b-spline (mazama.compgeom.b-spline/b-spline v+ v*
                                                      [[-4 -4] [-2 4] [2 -4] [4 4]]
                                                      knots)]
      (check-spline b-spline
                    knots
                    {0 [-4 -4]
                     1/2 [0 0]
                     1 [4 4]})))
  (testing "http://nurbscalculator.in/ 2"
    (let [knots [0 1 2 3 4 5 6 7]
          b-spline (mazama.compgeom.b-spline/b-spline v+ v*
                                                      [[-4 -4] [-2 4] [2 -4] [4 4]]
                                                      knots)]
      (check-spline b-spline
                    knots
                    {3 [-5/3 4/3]
                     7/2 [0 0]
                     4 [5/3 -4/3]})))
  (testing "http://nurbscalculator.in/ 3"
    (let [epsilon 1/1000
          control-points [[-4 -4] [-2 4] [2 -4] [4 4] [-0.7 -2.4] [1.3 -4.8]]
          knots [0 0 0 0.3 0.5 0.5 0.6 1 1 1]
          b-spline (mazama.compgeom.b-spline/b-spline v+ v* control-points knots)]
      (check-spline b-spline
                    knots
                    {0.3 [-0.88 -0.16]
                     0.5 [3.333 1.333]
                     0.6 [2.388 1.6]}
                    epsilon)))
  (testing "https://pages.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/de-Boor.html"
    (let [control-points [-1 1 3 4 6 9 0]
          knots [0 0 0 0 1/4 1/2 3/4 1 1 1 1]
          b-spline (mazama.compgeom.b-spline/b-spline + * control-points knots)]
      (check-spline b-spline
                    knots
                    {2/5 (shene-de-Boor + * control-points)}))))

;; https://www.cs.montana.edu/courses/spring2009/525/dslectures/CoxdeBoor.pdf
(deftest Starkey-test
  (testing "Example 1 - basis function"
    (let [knots [0 1 2 3 4 5]
          b-spline (mazama.compgeom.b-spline/b-spline v+ v*
                                                      [[1 1] [2 2] [3 0] [0 0]]
                                                      knots)]
      (check-spline b-spline
                    knots
                    {1 [1 1]
                     2 [2 2]
                     3 [3 0]
                     4 [0 0]})))
  (testing "Example 2 - basis function"
    (let [knots [0 0 1 2 3 3]
          b-spline (mazama.compgeom.b-spline/b-spline v+ v*
                                                      [[1 1] [2 2] [3 0] [0 0]]
                                                      knots)]
      (check-spline b-spline
                    knots
                      {0 [1 1]
                       1 [2 2]
                       2 [3 0]
                       3 [0 0]}))))
